function [kt1,kt]=Chapter4_1
% global vlast beta delta theta k0 kt
hold off
hold all
 vlast=zeros(1,100);
 k0=0.06:0.06:6;
 beta=.98;
 delta= .1;
 theta=.36;
 numits=240;
 for i=1:numits
     for j=1:100
         kt=j*.06;
         ktp1=fminbnd(@valfun,0.01,6.2);
         v(j)=-valfun(ktp1);
         kt1(j)=ktp1;
     end
     if i/48==round(i/48)
         plot(k0,v)
         drawnow
     end
     vlast=v;
 end
 hold off
 plot(k0,kt1)
 function val=valfun(k)
 g=interp1(k0,vlast,k,'linear');
 kk=kt^theta-k+(1-delta)*kt;
 if kk<=0
     val=-888-800*abs(kk);
 else
     val=log(kk)+beta*g;
 end
 val=-val;
 end
end